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ABSTRACT 

We stack WMAP 7-year temperature data around extragalactic point sources, showing 
that the p rofiles are consistent with WM AP's beam models, in disagreement with the 
findings of Sawangwit & Shanks ( 2010a ) . These results require that the source sample's 
selection is not biased by CMB fluctuations. We compare profiles from sources in the 
standard WMAP catalog, the WMAP catalog selected from a CMB-free combination 
of data, and the NVSS catalog, and quantify the agreement with fits to simple para- 
metric beam models. We estimate the biases in source profiles due to alignments with 
positive CMB fluctuations, finding them roughly consistent with those biases found 
with the WMAP standard catalog. Addressing those biases, we find source spectral 
indices significantly steeper than those used by WMAP, with strong evidence for spec- 
tral steepening above 61 GHz. Such changes modify the power spectrum correction 
required for unresolved point sources, and tend to weaken somewhat the evidence 
for deviation from a Harrison- Zel'dovich primordial spectrum, but more analysis is 
required. Finally, we discuss implications for current CMB experiments. 
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1 INTRODUCTION 

A telescope's beam or point spread function, which dampens 
the instrumental response to fluctuations on small scales, 
is prominent among the systematic effects that must be 
well-understood for Cosmic Microwave Background (CMB) 
measurements. The act of observing convolves the sky with 
the telescope beam, so observations of bright objects that 
are point-like (compared to the beam size) provide an 
obvious check on the beam pattern. The Wilkinson Mi- 
crowave Anisotropy Probe (WMAP) team used observations 
of Jupiter to measure the beams ( Jarosik et al. 20111 iHill 



et al.|2009||Jarosik et al.|2007||Page et al.|2003| ), but in prin- 
ciple any suitable objects will work, including extragalactic 
point sources. 



Sawangwit fc Shanks (2010a I used the 5-year WMAP 



data to construct stacked profiles around sources detected 



in the WMAP point source catalog (Wright et al. 20091 



Whitbourn et al. (20111 repeated this analysis with the 7- 



year data, using WMAP and external catalogs. Both papers 
report some intriguing discrepancies: for WMAP's differenc- 
ing assemblies at 40-90 GHz, they found substantial offsets 
at large scales, and the profiles appear broader than the 
beam patterns from the Jupiter measurements. They ad- 
dress and discard several possible explanations for the ef- 
fect: extended radio sources, source clustering, and selection 
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bias near the catalog threshold, and favor a nonlinearity in 
WMAP's response to Jupiter, which has a peak tempera- 
ture ~ 3 orders of magnitude higher than CMB fiuctuations 
at these frequencies. They note that a cosmological analysis 
based on the window functions computed from their stacked 
profiles (instead of the Jupiter model) would significantly 
change WMAP's basic cosmological results, for example by 
changing the height and location of the first acoustic peak in 



the power spectrum. In Sawangwit & Shanks (2010b I, the 



authors present this finding as a challenge to the ACDM 
paradigm. 

Despite these findings, several contrary lines of evidence 
indicate that the Jupiter-based beam models sufficiently rep- 
resent the true beam patterns for the WMAP telescope. 
First, WMAP and several pre- and post- WMAP experi- 
ments (from the ground and from balloons) have over the 
past decade found similar CMB power spectra on the scales 
where they overlap. Telescopes for which the beam scale 
differs substantially from WMAP are particularly power- 
ful probes of this consistency (see WMAP comparisons to 
QUaD : [Brown et al.l[2009| ACBAR: [Reichardt erari[2009l 



ACT: Hajian et al 



2010 and SPT: 



Keisler et al. 



20111 



Second, using WMAP data alone, comparison of the power 
spectra derived from the different differencing assemblies can 
be a useful test of errors in the beam models. A frequency- 
dependent beam systematic error of the size considered by 
[Sawangwit fc Shanks| ( [2010a[ l would almost certainly show 
up in estimates for the unresolved point source contribution 
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at high-L |Nolta et alT] ( |2009| ) and [HuflFenberger et al.| ( |2008 1 
extensively studied multifrequency combinations of Jupiter- 
beam-corrected power spectra from WMAP for this pur- 
pose, and found no such substantial beam anomaly. These 
methods constrain the relative beam window functions be- 
tween different assemblies at roughly the percent level. For 
example, a bump at Z < 200 in the WMAP 3-year data's un- 
resolved point source estimate disappeared with the beam 
revision of Hill et al. ( 2009[) for the 5-year data. This revi- 



sion changed the measurement of first acoustic peak's height 
by 2-3 percent, much smaller than the change in the beam 



proposed by Sawangwit & Shanks (2010a 


1 or Whitbourn 


et al. (20111. (See also Giannantonio et al. 


2010 for further 



discussion of the strengths of the ACDM model despite the 
misgivings of Sawangwit fc Shanks|2010b ) 

However, the issue with the stacked point sources 
remains: why should stacked extragalactic point sources 
present a profile so different from the telescope beam? This 
is a worthwhile question. Here, we explore it by stacking 



sources in the same way as Sawangwit & Shanks (2010a I, 



and we are able to reproduce their basic results for the 
WMAP standard catalog. However, by using alternative 
source catalogs, we by contrast find source profiles compat- 
ible with the WMAP beam models. We fit parameterized 
models to quantify the beam effects, and further explore sys- 
tematic selection biases. In section |2] we present our data 
selection and analysis methods, while in section[3]we discuss 
and interpret our results. Finally in section [4] we draw our 
conclusions. 



2 METHODS 
2.1 Data selection 

We base our analysis on the WMAP 7-year maps and point 
source catalogs, available from the LAMBDA website]^ In 
turn we use maps for the individual differencing assemblies 
(DAs), with and without a foreground template removed. 
We focus on Q-band (41 GHz), V-band (61 GHz), and 
W-band (94 GHz). All maps are at HEALPijj^ resolution 
A^sidc = 512 (6.9' pixels). 

We mask the sky to exclude pixels and sources near 
the galactic plane or other extended structures (like the 
LMC), while retaining bright sources away from extended 
foregrounds. We begin with the WMAP 7-yr temperature 
analysis mask. Then we invert the WMAP point source mask 
to retain those pixels near sources excluded by temperature 
analysis mask. The negative side effect is that some pixels 
are now kept around sources in high-foreground regions. To 
eliminate these, we smooth with a 2° FWHM Gaussian and 
apply a threshold at 90 percent of the smoothed maximum. 
This expands the mask slightly, eliminating the problem pix- 
els. Our resultant mask excludes the highest foreground re- 
gions and leaves 76.8 percent of the sky available for our 
analysis. 

The standard WMAP catalog contains 471 sources, but 
our mask excludes 38. Closely paired sources can spoil the 



profile, so we additionally require that the sources be iso- 
lated. Since we explore the source profiles out to Smax = 1.0°, 
we exclude all catalog sources which are separated by less 
than 20niax from another source. (Both members of the 
pair are excluded.) This prevents overlaps in computing our 
stacked profiles, and removes an additional 68 sources from 
our analysis (20 of these sources have another source even 
within 1.0°). After the two cuts, we have 365 sources on 
which we base our stacking analysis. The positional uncer- 
tainty in the W MAP catalog is 4' ( [Bennett et al.|2003[[CIien 
fc Wright|2008l ). 

The WMAP team also provides a catalog selected from 
CMB-subtracted maps using the multifrequency method of 



Chen & Wright (20081. The catalog contains 417 sources. 
Our mask eliminates 52 sources, then we cull 34 members 
of close pairs, leaving 331. The positional uncertainty in the 
CMB-free catalog is smaller, about 2'. Of these 331 sources, 
260 (or 78.5%) lie within 0.3° of a WMAP standard catalog 
source. 

Finally, we check our results by stacking on the source 



catalog from the NRAO VLA Sky Survey (NVSS, Condon 



et al. 1998 1, which surveyed the sky at 1.4 GHz for S > -40° . 
At that frequency, the CMB is not a significant component 
of the emission, and cannot affect the catalog selection. We 
cut at 2 Jy (at 1.4 GHz) to take the brightest 762 sources. 
Masking eliminates 350 sources, and culling close pairs elim- 
inates another 129 sources, leaving 283. Accounting for the 
reduced sky coverage, this gives a similar density of sources 
as the WMAP catalogs. The positional uncertainty of the 
bright NVSS sources is very small, < 1". Of these 283 NVSS 
sources, a total of 71 (or 25.1%) lie within 0.2° of a WMAP 
standard catalog source. Thus NVSS sources represent a sig- 
nificantly distinct population. 



2.2 Stacked profiles 

We stack our sources simply by defining angular annuli 
around reported catalog positions, then averaging the pixels 
whose centers fall into those annuli. 

The shape of the profile is subject to pixelization ef- 
fects, both from the angular binning into annuli and from 
the map pixelization. WMAP's beam models are tabulated 
at much higher angular resolution than the profiles, so must 
also be binned for direct comparison to the profile, or for 
the computation of ■ 

On an infinite resolution map, binning will simply av- 
erage the beam model (b) over annuli (indexed by i) to get 
the binned model (m): 



de e b{e) 



de 9. 



(1) 



^ http://lambda.gsfc.nasa.gov/ 



http : //healpix . j pi . nasa . gov 



The angular weighting causes the peak to be slightly sup- 
pressed. 

For point sources in realistic maps, the map pixel size 
and shape also influence the profile, broadening and blunt- 
ing the peak. This affects the smallest beams the most. To 
simulate this effect, we placed 1000 synthetic sources with 
the shape of each WMAP beam at random locations in a 
high-resolution HEALPix map, then reduced the resolution 
to A'aide = 512 and computed a stacked profile around the 
source center. (We used A'^sidc = 4096 for the high-resolution 




Angular separation (degrees) 



Figure 1. Smearing and pixelization effects compared to tlie 
input beam model for the W3 DA. After convolution with the 
HEALPix pixel window function, we show the binning of the 
input beam into annuli, where horizontal bars indicate the bin 
width. This is very close to the stacked profile from a Monte 
Carlo of synthetic sources. We also plot the beam model con- 
volved with a = 2' or 4' Gaussian to represent uncertainty in 
the source position. 
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Angular separation (degrees) 

Figure 2. The covariance matrix due to CMB and noise fluc- 
tuations for a profile in the W3 differencing assembly, using the 
standard WMAP catalog. Other DAs are qualitatively similar. 
CMB fiuctuations correlate across the range of angular separa- 
tions, while white noise contributes only to the diagonal, and is 
suppressed by the larger number of observations in annuli further 
from the center. 



map; using A^sidc = 2048 makes only a minor diflterence at 
W band, and is immaterial for V and larger beams.) This 
results in a profile which is notably blunted at the peak, 
and then larger than the input beam at 0.2° < < 0.4°. At 
the peak, pixel effects suppress the W band beams slightly 
more than ten percent, V band by about five percent, Q by 
three percent, and the K and Ka band profiles by less than 
two percent. Normalizing the pixelized, blunted profiles to 
the unbinned beam model at the peak can make the binned 
profile's tail look too heavy, but this alone is insufficient to 
account for the broad profiles seen by|Sawangwit fc Shanks| 
(2010a I. The source profiles in maps are well-reproduced by 



convolving the beam with the map pixel window function 
and then binning. This is the strategy we use later for our 
parametric beam models. We illustrate these binning and 
pixel effects in Figure[T]for the W3 DA. The W-band beams 
are WMAP's narrowest in the main lobe and this particular 
beam has an interesting shape, with a shoulder from 0.2- 
0.4 degrees. We model the uncertainty in source position by 
convolving the profile with a Gaussian of appropriate size 
(a = 4' for the standard WMAP catalog, a = 2' for the 
CMB-free catalog, and a = for the NVSS catalog). 

We consider two contributions to the covariance matrix 
for these profiles, from detector noise and from background 
CMB fluctuations. Our final covariance matrix is the sum of 
these components. Below we use these covariance matrices 
when minimizing for model fitting. We not make any cor- 
rection for the finite number of sources stacked in the profile, 
which will slightly modify the errors because of positional 
uncertainty (deviating from the Gaussian convolution kernel 
appropriate for an infinite number of sources). 

For non-overlapping source profiles, under the assump- 
tion of white noise, the detectors contribute a diagonal com- 
ponent to the covariance. This term we compute analytically 



from the noise variance per pixel provided by the WMAP 
team for each differencing assembly's map (Appendix [A|). 

Under the assumption that source positions are uncor- 
related with CMB fluctuations, we can construct the covari- 
ance due to the background CMB. For the WMAP standard 
catalog, we will later see that this assumption is unsound. 
The covariance can be computed analytically, but is more 
practically computed with a Monte Carlo method, as fol- 
lows. 

For the CMB power spectrum, we use the WMAP best- 
flt ACDM model. We compute stacked proflles on CMB- 
only simulations (including the beam and pixel window 
functions) and combine them to produce our estimate. We 
use 1600 Monte Carlo realizations, although using half that 
number changes our typical estimate of by just one half 
of one percent. 

Unlike the noise, the covariance due to CMB fluctua- 
tions is strongly correlated between proflle bins (Figure [2|. 

The covariance matrices for the standard WMAP cat- 
alog and the WMAP CMB-free catalog are very similar be- 
cause they share a large number of sources. The covariance 
for NVSS sources is somewhat larger for two main reasons. 
First, the number of NVSS sources is smaller. Second, NVSS 
did not survey the South ecliptic pole (5 ~ —66.6). WMAP 
noise is notably suppressed at the ecliptic poles due to a 
higher number of observations, so the NVSS proflles are 
taken from higher noise portions of the sky. This increases 
the errors accordingly below. 



2.3 Minimizing for Model Fitting 

Below we explore parametric models for the proflle, typically 
applying an amplitude factor and an offset. To flt a model 
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to our stacked source profile we minimize: 

(a) = ^ [P. - m, {a)]Cr^ ' [P, - m, {a}] . 



(2) 



Here P is the stacked profile and m is a binned model for 
the profile, which in turn is based on the WMAP beam, 
the HEALPix pixel window function, and the positional un- 
certainty appropriate for the catalog. The indices i,j run 
over the angular bins. This model additionally depends on 
a set of parameters a (like an amplitude, etc.). When the 
model depends on the parameters linearly, we can solve for 
the best fit parameters and their covariance matrix alge- 
braically. Otherwise we use an implementation of Powell's 
direction sets method to minimize for nonlinear models. 
We compute parameter covariances with 



[Cov(, 



ov Q!„, a, 



(3) 



although this is strictly only applicable when the fit is good. 



3 RESULTS 

3.1 Stacked source profiles 

For WMAP catalog sources, we show our stacked profiles 
from the DA maps in Figure [3] All plots are in thermody- 
namic temperature and we leave off the error bars to avoid 
crowding the plot, but show them below. Consistent with 
the spectral energy distribution typical of these sources, the 
average profile is brighter in Q band than in V band, and 
brighter in V than in W. Below we quantify that the pro- 
files are significantly wider than the Jupiter-modeled beams. 
The heavy tails in the profiles appear similar in V and W at 
d > 0.4°, as if the beam profiles are sitting atop a common 
fiuct nation. 

We see large scale offsets out to 1 degree, the same as 
Sawangwit fc Shanks] (|2010a[). The offset is 30-40 /xK in 



Q band and 10-20 fiK in V and W. ISawan gwit fc Shanks| 
( 2010a I attribute this offset to large scale CMB fiuctuations, 



and fit it at large angular scales. We disagree with this in- 
terpretation. CMB fluctuations would not cause the sub- 
stantially larger offset in Q band. Since the same sources 
are stacked in each band, the underlying CMB fiuctuations 
should be the same, up to the effects of beam smoothing. 
That Q band is so much higher suggests that galactic fore- 
grounds with a steeply falling spectrum, like synchrotron, 
might be responsible for this offset. Indeed the profiles com- 
puted in the foreground-reduced maps do not show the large 
scale offsets. Here, the Q band profiles too seem to be sitting 
atop the common fiuctuation, joining V and W at 9 > 0.5°. 

Despite some useful features, the foreground-reduced 
maps are inappropriate for probing the stacked profiles in 
detail. The templates for foreground removal also subtract 
away a portion of each source's flux. One foreground tem- 
plate is constructed from the K and Ka bands, and when the 
template is subtracted it creates a depression around each 
source which is the size of the larger beam scale in those two 
channels. 

The WMAP CMB-free catalog also selects sources that 
are bright in the WMAP bands; nearly 80 percent of the 
sources are the same. These stacks (Figure |4]| have similar 



peak temperatures to the WMAP sources at Q-band, and 
are slightly lower in V and W. Furthermore, the tails of 
these profiles are much less heavy than for sources from the 
standard catalog. 

The stacked profiles around NVSS catalog sources show 
that on average the NVSS sources are dimmer than the 
WMAP catalog sources, and the signal-to-noise is lower. 
This is consistent with their selection at lower frequency (at 
1.4 GHz), so that more falling spectrum sources and fewer 
GHz-peaked sources are represented here. As expected, the 
same pattern of source peak temperature dropping from Q 
to V to W holds. 



3.2 Parametric beam fits 

We fit a two-parameter model to the stacked profiles: 



m{e) = Sb{e) + mo 



(4) 



where S is the amplitude, b{9) is the Jupiter-modeled 
WMAP beam, normalized and smoothed with the pixel win- 
dow function and Gaussian positional uncertainty, and mo 
is an offset to represent large scale foreground or CMB con- 
tamination. Comparisons to the data use the model after 
binning. We enforce model beam normalizations so that the 
2D integrals under beams are unity. 



2n de e b{e) = i. 



(5) 



Thus b{9) has the units of inverse solid angle. Therefore the 
amplitude S has units of temperature times solid angle, or 
equivalently flux density. For each DA, conversions to fiux 
densities (dB/dT) are computed at the effective frequencies 
given in Page et al. ( 2003 1 for synchrotron-type falling spec- 



trum emission. 

We first examine the best-fit models for the W3 DA; 
the behavior is similar for the other DAs in the V and W 
bands. Figure [5] shows the best-fit model for this DA and 
the model's residuals for each of the three catalogs. 

The stack on the WMAP standard catalog shows a 
heavy tail compared to the fit, with large positive resid- 
uals. One seeming peculiarity of this fit is that the fitted 
model is below the data at most scales, and at least the 
model's peak could be brought closer to the data by in- 
creasing the amplitude. Such a change nonetheless worsens 
the of the fit, and the strong off-diagonal components 
in the covariance matrix caused by CMB fluctuations (Fig- 
ure [2]| are responsible. Because, unlike the stacked data, the 
WMAP beam models do not have such heavy tails, the flt 
will fall below the data at 0.5 degrees for any reasonable am- 
plitude. Different bins are so strongly correlated that this in 
turn causes the x^ fit to prefer that all bins be below the 
stacked profile, including those at ^ < 0.2°. Boosting the 
amplitude to raise the model near the peak, which naively 
would appear to improve the flt, actually worsens x^- Be- 
cause of the bin-to-bin correlations, in this case the fit 
prefers, in order: (1) all bins low, (2) only some bins low, 
(3) some bins high and some bins low. The other DAs show 
this same behavior stacked on the WMAP standard catalog, 
except for Q-band fltted to the stack from the foreground- 
reduced maps where the residuals are negative and caused 
by the foreground template over-subtraction (and the same 
reasoning regarding correlated bins applies). 



Stacking catalog sources in WMAP data 5 




0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 



Angular separation 9 (deg) Angular separation 9 (deg) 

Figure 3. Stacked profiles around sources from the WMAP point source catalog for all DAs, where the horizontal bars show the angular 
bin width. The source amplitude decreases from Q-band (red circles) to V-band (green triangles) to W-band (blue squares). The best-fit 
models for Ql, VI, and Wl are shown as dotted lines, and are poor fits to the profiles. Left: using the raw WMAP maps for each 
differencing assembly. Right: the same, but using foreground-reduced maps. 
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Figure 4. Similar to left panel of Figure |3] but stacking on the alternative source catalogs that are selected from maps without a 
significant CMB contribution. The beam models fit much better here. Left: WMAP CMB-free catalog. Right: NVSS catalog. 



By contrast, the profiles from the CMB-free and NVSS 
catalogs and do not show heavy tails. Remarkably, the shoul- 
der in the W3 beam appears to be recovered in the CMB- 
free profile. The NVSS source profiles have lower signal to 
noise, but are similar. The other DAs show the same for the 
CMB-free and NVSS catalogs. 

Thus, sources selected from maps that contain CMB 
produce stacked profiles with a bias, while sources selected 
from maps without a significant CMB contribution do not. 



These results differ from what Sawangwit & Shanks (2010a I 



and Whitbourn et al. (20111 report for their CMB-free and 



NVSS stacks. Below we find this bias is roughly consistent 
with the expectations of source selection bias due to CMB 
fiuctuations. 

This bias for the WMAP standard catalog causes the 
of the fit to be very poor in every DA, seen in Table [l] 
which also shows the fitted parameters and the probability 



to exceed by chance for 18 degrees of freedom (20 bins 
with two parameters). 



The CMB-free catalog probabilities are much more rea- 
sonable, except at Q band, where they are a little low. The 
probabilities for the NVSS catalog are everywhere reason- 
able. 



Despite some minor complications which we discuss be- 
low, we find no compelling evidence from this modeling that 
the the Jupiter-based WMAP beams are radically insuffi- 
cient to explain the source profiles, and we conclude that 
WMAP's beams for V- and W-bands, which are used for 
the cosmological analysis, are sound at this level. Further- 
more, the discrepancy with the profiles from the standard 
WMAP catalog bears the mark of source selection bias due 
to CMB fluctuations. We explore this below. 
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Figure 5. Stacked profile of catalog sources on the W3 WMAP map, fit with the amplitude-offset model (which includes the pixel 
window function and positional uncertainty). Plotted errors are from the diagonal of the covariance matrix. Top left: WMAP standard 
catalog sources, where the profile shows a bias. Top middle: WMAP CMB-free catalog sources. Top right: NVSS catalog sources. Bottom: 
binned residuals for each fit. 
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Figure 6. Like Figure|5] stacked profile of catalog sources but on the Ql WMAP map. Top left: WMAP standard catalog sources, where 
the profile shows a bias. Top middle: WMAP CMB-free catalog sources, where = 35.41 for 18 degrees of freedom is particularly high 
for this catalog. Top right: NVSS catalog sources. Bottom: binned residuals for each fit. 



3.3 CMB-free catalog residuals 

The CMB-free catalog profiles have some features that war- 
rant further discussion. For example, for every CMB-free 
profile, the first bin is low compared to the models for each 
DA (Figure [3] left panel). 

Figure |6] depicts the profiles and fits for Ql, which is 
the poorest fit for the CMB-free catalog. Here the residuals 
show some interesting patterns. In both Q-band DAs, the 
model for the CMB-free catalog exceeds the data slightly 
at 0.3° < 6 < 0.5°, leading to negative residuals. One 
possibility for this relates to the positional uncertainty. 
Brighter sources will have more accurate measured posi- 
tions, so stacked profiles, weighted toward the bright sources, 
should have less effective positional uncertainty than the cat- 
alog overall (quoted as 2' for this case). In addition, the mix 
of bright sources in each stack changes from band to band 
based on the frequency dependence of each source, so even 



considering a single catalog, the effective positional uncer- 
tainty can be different between the bands. 

When we allow the positional uncertainty to float in the 
Q-band fit, the minimum has positional uncertainty less 
than 2', but it is not consistent between Ql, where the fit 
prefers zero positional uncertainty, and Q2, which prefers 
1.9', and the probability to exceed x^ is still low, never ex- 
ceeding a couple percent (now for 17 degrees of freedom). 
Changing the positional uncertainty from 0' to 2.5' leads to 
a 5 percent increase in the inferred source amplitude. Fur- 
ther letting the model adjust the positional uncertainty in 
a flux-dependent way makes the model so flexible that it is 
difficult to draw any conclusions. 

The other peculiar feature in Figure[6]for the CMB-free 
Ql residuals is the alternating low- high pattern in the first 
several bins, and the first bin is quite low. This is visible 
in the first five bins in each DA except for Wl and W4, 
where only the first three and two bins (respectively) follow 
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WMAP catalog 



DA 


Amp. (Jy) 


Offset (fiK) 


v2 

A^=18 


Pi> x') 




Ql 


1.51 ± 0.01 


32.80 ± 4.07 


57.08 


6.01 X 10- 


6 


Q2 


1.51 ± 0.01 


33.45 ± 4.17 


44.24 


5.34 X 10" 


4 


VI 


1.46 ± 0.02 


12.49 ± 4.27 


49.45 


9.12 X 10" 


5 


V2 


1.43 ± 0.02 


13.97 ± 4.14 


67.23 


1.32 X 10" 


7 


Wl 


1.38 ± 0.04 


16.57 ± 4.33 


56.09 


8.63 X 10- 


6 


W2 


1.43 ± 0.05 


16.92 ± 4.27 


69.83 


4.83 X 10" 


8 


W3 


1.41 ± 0.05 


16.39 ± 4.19 


58.66 


3.36 X 10- 


6 


W4 


1.38 ± 0.04 


16.54 ± 4.21 


47.35 


1.88 X 10- 


4 


WMAP CMB-free catalog 


DA 


Amp. (Jy) 


Offset inK) 




P(> x') 




Ql 


1.50 ± 0.01 


36.97 ± 4.28 


35.41 


8.39 X 10- 


3 


Q2 


1.49 ± 0.01 


36.58 ± 4.27 


30.34 


0.034 




VI 


1.31 ± 0.02 


15.97 ± 4.26 


25.52 


0.111 




V2 


1.27 ± 0.02 


15.56 ± 4.39 


19.22 


0.378 




Wl 


1.04 ± 0.03 


17.32 ± 4.27 


16.86 


0.533 




W2 


1.05 ± 0.04 


16.85 ± 4.40 


24.60 


0.136 




W3 


1.03 ± 0.04 


16.25 ± 4.35 


17.16 


0.512 




W4 


1.05 ± 0.04 


17.23 ± 4.32 


8.89 


0.962 




NVSS catalog 


DA 


Amp. (Jy) 


Offset inK) 


A^=18 


P(> x') 




Ql 


0.77 ± 0.01 


31.04 ± 6.72 


17.05 


0.520 




Q2 


0.74 ± 0.01 


28.86 ± 6.72 


24.59 


0.137 




VI 


0.63 ± 0.02 


10.27 ± 6.57 


14.40 


0.703 




V2 


0.66 ± 0.02 


12.14 ± 6.73 


12.43 


0.825 




Wl 


0.58 ± 0.03 


13.17 ± 6.70 


12.02 


0.846 




W2 


0.54 ± 0.04 


12.81 ± 6.57 


15.92 


0.598 




W3 


0.59 ± 0.04 


14.04 ± 6.77 


24.41 


0.142 




W4 


0.54 ± 0.04 


13.81 ± 6.71 


9.25 


0.954 





Table 1. Parameter values and the goodness-of-fit for a simple 
amplitude-offset model based on the WMAP beam, for stacked 
profiles from the three different catalogs. The probability to ex- 
ceed indicates that the model fits the WMAP catalog profiles 
poorly, the CMB-free catalog profiles much better (although Q 
band's probability is low), and the NVSS profiles well. 



the pattern. This cannot be ringing in the maps due to the 
source signal, because the pattern is absent in the NVSS 
profiles. Positional uncertainty or systematic mis-centering 
of the sources do not cause an alternating pattern. 

CMB fluctuations, which are common to all the DA 
maps (up to beam smoothing), are another candidate. The 
covariance matrix (Figure [2| shows that bins 1-3-5 are more 
strongly correlated to each other (30-35 fJ-K^) than they are 
to bins 2-4 (20-25 IJ-K^), so CMB fluctuations could conceiv- 
able have an effect shaped like this (with the offset removing 
the fluctuation common to all bins) and be common across 
the bands. However, such an effect should only appear at 
the couple of fjK level, smaller than what is seen here, and 
would not be particularly stronger in Q. 

Slight underestimates of the covariance due to posi- 
tional uncertainty and the finite number of sources, which 
would be stronger in Q, are another possibility, but we have 
not explored it in detail. 



3.4 Flux density and spectral index 

The source amplitudes are similar in Q band for the WMAP 
standard and CMB-free catalogs in Table [l] but the ampli- 



Catalog 


"QV 




"QW 


WMAP (*) 


-0.11 ± 0.03 


-0.07 ± 0.04 


-0.09 ± 0.02 


CMB-free 


-0.36 ± 0.03 


-0.50 ± 0.05 


-0.43 ± 0.02 


NVSS 


-0.40 ± 0.06 


-0.31 ± 0.09 


-0.35 ± 0.04 



Table 2. Frequency scaling for the three catalogs. The errors rep- 
resent the error on the catalog averages propagated from Table [T] 
and not the much larger intrinsic scatter for sources. 
(*) The WMAP standard catalog values are biased by the poor 
fits and shown only for comparison. 



tudes in V and W are substantially less in the CMB-free 
catalog. The selection bias causes the sources in V and W 
to appear much too bright, and this has a significant effect 
on the measured frequency dependence of sources. 

From the best-fitted source amplitudes, we can examine 
the scaling for the mean source flux density, writing 

Sv/Sq = (i^v/j^q)"«^ (6) 

and so forth. These are displayed in Table [2] For the WMAP 
standard catalog, where we know the flux densities to be bi- 
ased, we find a mean spectral index a ~ —0.09, the same 



as Wright et al. (20091 found on a source-by-source basis. 
However, the spectral index from the CMB-free catalog is 
much steeper and shows steepening above 61 GHz, from 
aqv ~ —0.36 ±0.03 to avw ~ —0.50 ±0.05, a difference sig- 
nificant at the ~ 5a level. The NVSS sources have similarly 
steeper indices than the WMAP standard catalog sources, 
but there is no significant change in the spectral index with 
frequency. 

These frequency scalings are employed in the estimation 
of the power spectrum of unresolved point sources, which is 
necessary to correct the CMB power spectrum. The source 
correction effectively changes the tilt of the measured power 
spectrum. WMAP's CMB power spectrum results, which 
are the basis of the cosmological measurement, employ a 
combination of V- and W-band data only. 

The point source correction uses multifrequency esti- 
mators that scale the unresolved sources from Q to the V 
and W bands. If the power in Q is held constant, steep- 
ening the spectral index according to the CMB-free cata- 
log values reduces the required correction in V by about 
20 percent and in W by about 40 percent. Comparing to 



Huffenberger et al. ( 2006 2008 1 such changes should raise 



the value of the scalar index for primordial perturbations 
(ris) by ~ 0.01 or so, which would affect the statistical sig- 
nificance for WMAP's confirmation of the inflationary pre- 



diction that Us < 1 (see Komatsu et al. 2011 which gives 
Us = 0.968 ± 0.012). However, some caution is appropriate: 
the unresolved source contamination is slightly smaller in 
the more recent WMAP release, making it less susceptible, 
and this is a substantial enough change to the scaling that a 
more complete analysis is really required to be quantitative. 

3.5 Selection bias from CMB anisotropy 

The same fluctuations that are responsible for biases in the 



catalog selection and source counts (Eddington 19131 will 
bias the source proflles. For the WMAP catalogs, the sources 
are mostly selected at lower frequency, where they tend to 
be brightest, and there is a signiflcant mismatch between 
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the beam scales at K-band (49' FHWM) and th e Q, V, and 
W bands (29', 20', 13' FWHM, respectively) ([Page et al. 



2003). Although the noise is distinct, the same CMB is seen 
by each of the DAs (up to beam smoothing). The profiles 
in the higher-frequency DAs can therefore be broadened by 
CMB fluctuations that cannot be distinguished from source 
flux in the lower frequency maps. 

To roughly quantify the bias as a function of source flux 
density, we performed Monte Carlo simulations, injecting 
sources of known flux into synthetic realizations of the CMB 
and detector noise and then finding these sources. 

FollowinR the WMAP source detection procedure 
(Wright et al. 20091, the maps are weighted by (A'oba)^''^, 



then filtered in harmonic space by bi/{bi'wfC[ + Ni) to 
maximize the signal-to-noise of the recovered sources. The 
noise spectrum is computed from the pixel noise, pixel area, 
and the mean over the sky of the inverse number of observa- 
tions: Ni = (Tof2pixAf~j^^. The pixel noise factors ag are given 
by the WMAP data release for each DA. 

The synthetic sources that exceed 5cr in the filtered 
map are taken to exceed the catalog threshold. Our syn- 
thetic source selection is slightly simpler than the WMAP 
selection. First, we select only in K (where most sources are 
brightest), instead of the five-band selection used by WMAP. 
Second, we use the center of the local maximum pixel as 
the source position, rather than fitting for the position of 
each source candidate. We stack the unfiltered maps around 
the recovered sources and examine the source profiles. We 
show the bias for several input source fiux densities in Fig- 
urejT] For a 1.0 Jy source the bias is a little less than 20 /iK 
at the source position and falls to half-maximum at about 
0.5 degrees. For faint sources, only those that fall on back- 
ground peaks are recovered, so the bias is large. Almost all 
bright sources are found, but slight biases remain because 
of errors in the source positions, which are shifted to favor 
background peaks. 

After evaluation at several flux densities, we can inter- 
polate to approximate Pbias(^, S), the bias in the profile as a 
function of fiux density. Then the expected bias in a catalog 
can be estimated with 



-Pbias(^) — 



ds^p^Ue,s), (7) 



where dNcnt/dS gives the distribution of fiuxes in the cata- 
log, accounting both for the intrinsic source counts and the 
selection function, and A'cat is the total number of sources. 
We estimate this integral by evaluating a sum over sources 
in the WMAP CMB-free catalog, which is free from CMB 
contamination and lists noise-bias corrected source fluxes: 



Pbia.(^) 



iVca 



^ ^ -Pbias(^, Si 



(8) 



The median K-band flux in this catalog is 0.9 Jy and the 
minimum de-boosted flux is 0.1 Jy. The resulting expected 
bias is about 25 /iK at the peak and quite broad. Despite the 
crudeness of our estimate, which is slightly larger than the 
residual for the WMAP standard catalog profile in Figure [5j 
where the fit takes up some of the residual, the bias due to 
source selection appears to be the probable explanation for 
the broad profiles found here and in |Sawangwit fc Shanks| 
(2010a|l for the WMAP standard catalog. 



M 40 



0.5Jy 



0-7Jy 




).5 1.0 1.5 2.0 2.5 3.0 3.5 4.0 

Flux density S (Jy) 




0.2 0.4 0.6 

Angular separation (degrees) 



Figure 7. Profiles of CMB fluctuations at the location of found 
sources, for sources witfi input S = 0.5 Jy— 4.0 Jy in WMAP's Kl 
band. Tfiis represents a bias in the source profile due to the source 
selection. We use the same set of CMB and noise realizations in 
each case, accounting for some common fluctuations in the curves. 
Although ^ 100 percent of sources with S > 4.0 Jy are detected, 
slight biases remain due to position errors. The mean profile bias 
(Pljias(9)) is roughly what we expect for the set of sources in the 
WMAP catalog. 



4 CONCLUSIONS 

We stacked point sources from three different catalogs on the 
maps from eight of WMAP's differencing assemblies in the 
Q, V, and W bands. For the WMAP standard catalog, we 
see evidence for residual CMB fiuctuations that bias the pro- 
files. Some complications remain in the profiles for CMB-free 
catalog sources, but for V and W band the fits are reason- 
able. For NVSS sources, the fits are reasonable for all DAs. 
Therefore, when sources are selected from data that contain 
no significant CMB contribution, we find no compelling ev- 
idence that the beams for WMAP differ substantially from 
the Jupiter-based models. These conclusions directly contra- 



dict those of Sawangwit & Shanks (2010a I and Whitbourn 
|et al.| ( |2011| ), although they do not report the statistical 
significance of their result. The reason for the discrepancy, 
especially for the CMB-free and NVSS catalogs, is not clear. 
One possibility relates to the culling of close pairs, where we 



cut at 2 degrees, while [Whitbourn et al.| ( |2011[ ) cuts NVSS 
pairs at 1 degree and WMAP CMB-free pairs not at all. 
At the same time, they follow the profile out further than 
we do, beyond two degrees. Their method for background 
subtraction also differs from our offset fit. 

The biases in the profiles erroneously boost the inferred 
source fiux, especially for the smaller beams, and this affects 
the source spectral index. Using fiuxes from the CMB-free 
catalog, the spectral indices are significantly steeper, and 
show spectral steepening at high frequency. This in turn 
affects the point source subtraction for the power spectrum. 

The issues we confronted in this work (pixel and bin- 
ning effects, positional uncertainty, foreground subtraction, 
selection bias) are common to all microwave experiments 
that construct point source catalogs. Recent catalogs from 
ACT ( [Marriage et~aI1[20TT| ) and SPT ( [Vieira et ai:i[20Tot 
are much less susceptible to the CMB fluctuations because 
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the power in the CMB faUs so rapidly on arcminute angular 
scales. These catalogs take care to de-boost their measured 
flux densities for the biases due to noise and other sources. 
However, at these frequencies and angular scales, the back- 
ground of dusty galaxies can play a role similar to that the 
CMB plays here. Subsequent analysis of these catalogs could 
run into similar biases, for example, when stacking 220 GHz 
maps on 150 GHz source positions or vice versa, or stacking 
ACT 220 GHz sources in the SPT data where the beam is 
smaller. 



Planck data (Planck Collaboration et al. 2011al, with 
larger beams (5'-30 ), is more susceptible to CMB fluctu- 
ations, so the same cautions apply as for WMAP, both for 



the recent Early Source Catalog ( Planck Collaboration et al 



2011b I and for the flnal band-merged catalog. 



ACKNOWLEDGMENTS 

We thank Xi Chen, Reijo Keskitalo, Mike Nolta, and Ned 
Wright for useful comments and discussion. Some of the re- 
sults in this paper have been derived using the HEALPix 



package ( Gorski et al. 2005 1 and Healpy. This work was 



partially supported by a Space Education and Training 
Program award from the Florida Space Grant Consortium. 
KWS acknowledges the support of the Beyond the Book 
program from the University of Miami College of Arts and 
Sciences. KMH receives support from NASA-JPL subcon- 
tract 1363745. 



APPENDIX A: COVARIANCE FOR STACKED 
PROFILES 

The stacking procedure is equivalent to the linear operation: 



(Al) 



where Pi is the stacked profile in angular bin i, Mq is the 
gth pixel's value in the map, and a stacking operator Siq is 



Siq 



1 if pixel q is in anrmlus i about a source, 
otherwise. 



(A2) 

The catalog point source positions are implicit in the S ma- 
trix, thus the contraction of the map with the stacking oper- 
ator separately accumulates the contribution for each profile 
bin. 

If the map M is zero mean, as for CMB and detector 
noise, then the covariance in of profile bins is a simple func- 
tion of the covariance of maps. 



Cov(Pi,Pj) = 



(A3) 



For a map with only white detector noise, Uq, we have 
{uqUr) = 5qrao/Nohs,q- Thcrcfore, 

Cov(P„P,) = 



{HqS^q) 



2 '^q ^^Q^JI^obs 



obs,g 



(A4) 



If additionally we require that the sources are well- 
separated, that is, that the annuli of different sources never 
overlap (and the annuli around single sources naturally 
never overlap), then the bins are uncorrelated and the co- 
variance simplifies: 



Cov(Pi, Pj) = Sijao 



Eq ^i-l^iQ-^obs,, 



^Eq 'S'lq^ (Er ^^r) 

Since Siq has value zero or one, then SiqSiq — Siq, so 



Cov(Pi, Pj) = Sijao 



2 Eg Sig^obs 



(e,5,,) 



2 



(A5) 



(A6) 



Cij = Cov{Pi,Pj) = Sijal 



Eg ^iq^ohs.q 



Eg ^ig 



(noise only) 



(A7) 

where uoN^y^^ ^ is the noise variance in pixel q, provided 
by the WMAP team for each differencing assembly's map. 



Noting the similarity to the stacking operation (Eq. All, 
we compute this quantity using a slight modification to the 
stacking code, but applied to the A^J^^ map with a different 
normalization, and verify with Monte Carlo noise simula- 
tions. 



For the CMB, the term in angle braces in (A3 1 is simply 



the angular correlation function, obtained from the theory 
angular power spectrum Ci by a sum over Legendre polyno- 
mials: 

(MqMr) = Oj{eqr) = -!- V (2/ + 1) bf wf C I Pi (cOS 6 qr) . 

47r ^ — ' 

(A8) 

The beam window function for the map is represented by 
and the pixel window function by wf. In practice, however, 
it is computationally inconvenient to keep track of all possi- 
ble pixel separations (Oqr), and more convenient to compute 
this portion of the covariance using a Monte Carlo of CMB 
realizations. 
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